function [fun, grad] = PL_penalty_Nambu(beta, varargin)
% [fun, grad] = PL_penalty_Nambu(beta, varargin)
% 
% PL_penalty_Nambu computes the Nambu functional [1].
%
% See also: PL_penalty_parallel, PL_penalty_CTV, PL_penalty
%
% [1] N. Sochen, R. Kimmel, and R. Malladi, 
%     A General Framework for Low Level Vision,
%     IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 310–318, 1998
%
% 2013-07-25 --------------------------------------------------------------
% 
% Matthias J. Ehrhardt
% CMIC, University College London, UK
% 
% matthias.ehrhardt.11@ucl.ac.uk
% 
% -------------------------------------------------------------------------
    
    if iscell(varargin{1})
        K = length(varargin{1});
        input = varargin{1};
    else
        K = length(varargin);
        input = varargin;
    end
    
    % cell array for the gradients of all images
    G = cell(K,1);
    % cell array for the pairwise scalar products of the gradients
    S = cell(K, K);
    % cell array for the norms of the gradients
    N = cell(K, 1);   
    
    % For all pairs of images compute their scalar product of the
    % gradients. As the scalar product is symmetric it is sufficient to
    % skip almost half of the pairs.   
    for u = 1 : K        
        % calulate the gradient
        G{u} = 1./beta * Gradient(input{u}./beta, 'symmetric', 'central', 1, 0);
        % and store its norm
        N{u} = Norm_3(G{u}, 0);
    end
    
    for u = 1 : K
        for v = 1 : K
            if u ~= v
                if u < v
                    S{u, v} = Inner_prod_3(G{u},G{v});
                else
                    S{u, v} = S{v, u};
                end
            end
        end;
    end
                    
    c = numel(input{1});
    
    if nargout > 1
        grad = zeros(c* K,1);
    end
    
    eta = ones(size(input{1}));
    for u = 1 : K
        eta = eta + N{u}.^2;
        for v = 1 : K
            if v ~= u            
                eta = eta + 1./2 * (N{u}.^2 .* N{v}.^2 - S{u,v}.^2);
            end
        end
    end
    eta = sqrt(eta);

    if nargout > 1        
        for u=1 : K            
            I = (u-1) * c + (1 : c);
            kappa = ones(size(input{1}));
            for v = 1 : K
                if v ~= u
                    kappa = kappa + N{v}.^2;
                    tau = - S{u,v} ./ eta;
                    grad(I) = grad(I) - Diffusion(tau, input{v});
                end
            end
            
            kappa = kappa ./ eta;
            grad(I) = grad(I) - Diffusion(kappa, input{u});
        end
    end

    fun = sum(eta(:));
end